packages <- c("CIMseq", "CIMseq.testing", "tidyverse", "ggthemes")
purrr::walk(packages, library, character.only = TRUE)
rm(packages)

##DATA
load('../data/CIMseqData.rda')
load('../data/sObj.rda')
if(!dir.exists('../figures')) dir.create('../figures')
p <- plotUnsupervisedClass(cObjSng, cObjMul, palette('total'))
p

ggsave(
  plot = p,
  filename = '../figures/MGA.enge.classes.pdf',
  device = cairo_pdf,
  height = 180,
  width = 180,
  units = "mm"
)
p <- plotUnsupervisedMarkers(
  cObjSng, cObjMul,
  c("Lgr5", "Ptprc", "Chga", "Dclk1", "Alpi", "Slc26a3", "Atoh1", "Lyz1"),
  pal = RColorBrewer::brewer.pal(8, "Set1")
)
ggsave(
  plot = p,
  filename = '../figures/MGA.enge.markers.pdf',
  device = cairo_pdf,
  height = 180,
  width = 180,
  units = "mm"
)

p <- plotUnsupervisedMarkers(cObjSng, cObjMul, "Mki67")
ggsave(
  plot = p,
  filename = '../figures/MGA.enge.Mki67.pdf',
  device = cairo_pdf,
  height = 180,
  width = 180,
  units = "mm"
)

p <- plotUnsupervisedMarkers(cObjSng, cObjMul, "Hoxb13")
ggsave(
  plot = p,
  filename = '../figures/MGA.enge.Hoxb13.pdf',
  device = cairo_pdf,
  height = 180,
  width = 180,
  units = "mm"
)
classHeatmap(
  data = data.frame(
    gene = markers$gene,
    class = markers$cluster
  ),
  counts.log = getData(cObjSng, "counts.log"),
  classes = tibble(
    sample = colnames(getData(cObjSng, "counts")), 
    class = getData(cObjSng, "classification")
  )
)

Cell type markers

markers <- c(
  "Lgr5", "Ptprc", "Chga", "Dclk1", "Alpi", "Slc26a3", "Atoh1", "Lyz1", "Mki67",
  "Hoxb13", "Plet1", "Osr2", "Reg4", "Sval1"
)

gene.order <- c(
  "Hoxb13", "Osr2", "Atoh1", "Reg4", "Sval1", "Plet1", "Mki67", "Lgr5", "Alpi", 
  "Slc26a3", "Lyz1", "Dclk1", "Chga", "Ptprc"
)

d <- getData(cObjSng, "counts.cpm")[markers, ] %>%
  t() %>%
  as.data.frame() %>%
  rownames_to_column("Sample") %>%
  as_tibble() %>% 
  gather(gene, cpm, -Sample) %>%
  inner_join(tibble(
    Sample = colnames(getData(cObjSng, "counts")), 
    Classification = getData(cObjSng, "classification")
  )) %>%
  mutate(Classification = parse_factor(Classification, levels = classOrder.MGA('total'))) %>%
  mutate(gene = parse_factor(gene, levels = gene.order)) %>%
  arrange(Classification, gene) %>%
  mutate(Sample = parse_factor(Sample, levels = unique(Sample))) %>% 
  mutate(type = case_when(
    str_detect(Classification, "^S") ~ "Small intestine",
    str_detect(Classification, "^C") ~ "Colon",
    TRUE ~ "Other"
  )) %>%
  mutate(type = parse_factor(
    type,
    levels = c("Small intestine", "Colon", "Other")
  ))
## Joining, by = "Sample"
under <- d %>%
  mutate(idx = 1:nrow(.)) %>% 
  group_by(Classification) %>% 
  summarize(
    min = (min(idx)), max = (max(idx) - 1), median = median(idx)
  ) %>%
  ggplot() +
  geom_segment(
    aes(x = min, xend = max, y = 0, yend = 0, colour = Classification)
  ) +
  geom_text(
    aes(median, y = 0, label = Classification),
    angle = 90, nudge_y = -0.01, hjust = 1, colour = "grey10"
  ) +
  scale_colour_manual(values = palette('total')[classOrder.MGA('total')]) +
  ylim(c(-1, 0)) +
  theme_few() +
  scale_x_continuous(expand = c(0, 0)) +
  theme(
    axis.text = element_blank(),
    axis.title = element_blank(),
    axis.ticks = element_blank(),
    panel.border = element_blank()
  ) +
  guides(colour = FALSE)

over <- d %>% 
  ggplot() + 
  geom_bar(aes(Sample, cpm, fill = Classification), stat = "identity") + 
  facet_grid(gene ~ type, scales = "free", space = "free_x") +
  scale_fill_manual(values = palette('total')[classOrder.MGA('total')]) +
  theme_few() +
  theme(
    axis.text.x = element_blank(),
    axis.title.x = element_blank(),
    axis.ticks.x = element_blank(),
    strip.background.x = element_rect(fill = "white")
  ) +
  labs(y = "CPM") +
  guides(fill = FALSE)

library(patchwork)

p <- over + under + plot_layout(ncol = 1, heights = c(2, 1))
p

ggsave(
  plot = p,
  filename = '../figures/MGA.enge.geneBar.pdf',
  device = cairo_pdf,
  height = 300,
  width = 300,
  units = "mm"
)

Stemness dotplot

load('../data/seurat.rda')
markers <- c("Wnt3", "Wnt2b", "Wnt4", "Wnt5a", "Lgr5", "Axin2", "Ascl2")
gene.order <- markers

#scale.func <- scale_size
scale.func <- scale_radius
scale.min <- NA
scale.max <- NA

p <- getData(cObjSng, "counts.cpm")[markers, ] %>%
  t() %>%
  as.data.frame() %>%
  rownames_to_column("Sample") %>%
  as_tibble() %>% 
  gather(gene, cpm, -Sample) %>%
  inner_join(tibble(
    Sample = colnames(getData(cObjSng, "counts")), 
    Classification = getData(cObjSng, "classification")
  )) %>%
  mutate(type = case_when(
    str_detect(Classification, "^S") ~ "Small intestine",
    str_detect(Classification, "^C") ~ "Colon",
    TRUE ~ "Other"
  )) %>%
  group_by(gene, Classification) %>%
  summarize(mean = mean(cpm), pct = 100 * (length(which(cpm != 0)) / n())) %>%
  mutate(scaled.mean.exp = scale(mean)) %>%
  ungroup() %>%
  mutate(type = case_when(
    str_detect(gene, "^Wnt") ~ "Wnt ligand",
    gene %in% c("Lgr5", "Axin2", "Ascl2") ~ "Wnt target",
    TRUE ~ "error"
  )) %>%
  mutate(Classification = parse_factor(Classification, levels = classOrder.MGA('total'))) %>%
  mutate(gene = parse_factor(gene, levels = rev(gene.order))) %>%
  mutate(type = parse_factor(type, levels = c("Wnt ligand", "Receptor", "Wnt target"))) %>%
  ggplot() +
  geom_point(
    aes(Classification, gene, size = pct, colour = scaled.mean.exp)
  ) +
  facet_grid(type ~ ., space = "free", scales = "free") +
  scale_color_gradient(low = "white", high = "darkred") +
  scale.func(range = c(0, 18), limits = c(scale.min, scale.max)) +
  guides(
    size = guide_legend(
      title = "% expressed", title.position = "top", title.hjust = 0.5
    ),
    colour = guide_colorbar(
      title = "Scaled mean expression", title.position = "top", 
      title.hjust = 0.5, barwidth = 10
    )
  ) +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1), 
    axis.title = element_blank(), 
    legend.position = "top",
    legend.justification = "center"
  )
## Joining, by = "Sample"
p

ggsave(
  plot = p,
  filename = '../figures/MGA.enge.stemness.dotplot.pdf',
  device = cairo_pdf,
  height = 160,
  width = 250,
  units = "mm"
)

Relative frequency of cell types in singlets vs. deconvoluted multiplets

singlets <- c(table(getData(cObjSng, "classification")))
singlets <- singlets / sum(singlets)
deconv <- colSums(adjustFractions(cObjSng, cObjMul, sObj))
deconv <- deconv[match(names(singlets), names(deconv))]
deconv <- deconv / sum(deconv)
if(!identical(names(singlets), names(deconv))) stop("name mismatch")

p <- tibble(
  class = names(singlets),
  singlet.freq = singlets,
  multiplet.freq = deconv
) %>%
  ggplot() +
  geom_point(aes(singlet.freq, multiplet.freq, colour = class), size = 3) +
  scale_colour_manual(values = palette('total')[order(names(palette('total')))]) +
  xlim(min(c(deconv, singlets)), max(c(deconv, singlets))) +
  ylim(min(c(deconv, singlets)), max(c(deconv, singlets))) +
  geom_abline(slope = 1, intercept = 0, lty = 3, colour = "grey") +
  labs(x = "Singlet relative frequency", y = "Multiplet relative frequency") +
  guides(colour = guide_legend(title = "Cell Type")) +
  theme_few()

p

ggsave(
  plot = p,
  filename = '../figures/MGA.enge.sngMulRelFreq.pdf',
  device = cairo_pdf,
  height = 160,
  width = 250,
  units = "mm"
)
plotSwarmCircos(
  sObj, cObjSng, cObjMul, weightCut = 0, classOrder = classOrder.MGA('total'), 
  classColour = palette('total')[classOrder.MGA('total')], h.ratio = 0.85
)
## Joining, by = "class"

Detected duplicates - quadruplicates.
Only ERCC estimated cell number max 4.
Weight cutoff = 10.

# adj <- adjustFractions(cObjSng, cObjMul, sObj, binary = TRUE)
# samples <- rownames(adj)
# rs <- rowSums(adj)
# keep <- rs == 2 | rs == 3 | rs == 4

plotSwarmCircos(
  sObj, cObjSng, cObjMul, weightCut = 10, 
  classOrder = classOrder.MGA('total'), classColour = palette('total')[classOrder.MGA('total')], h.ratio = 0.85,
  alpha = 1e-3
)
## Joining, by = "class"

pdf('../figures/MGA.enge.circos.pdf', width = 12, height = 10, onefile = FALSE)
plotSwarmCircos(
  sObj, cObjSng, cObjMul, weightCut = 10, 
  classOrder = classOrder.MGA('total'), classColour = palette('total')[classOrder.MGA('total')], h.ratio = 0.85,
  alpha = 1e-3
)
## Joining, by = "class"
invisible(dev.off())